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We have proposed a method for the dynamic simulation of a collection of self-propelled 
particles in a viscous Newtonian fluid. We restrict attention to particles whose size and 
velocity are small enough that the fluid motion is in the creeping flow regime. We have 
proposed a simple model for a self-propelled particle, and extended the Stokesian Dy- 
namics method to conduct dynamic simulations of a collection of such particles. In our 
description, each particle is treated as a sphere with an orientation vector p, whose loco- 
motion is driven by the action of a force dipole S p of constant magnitude So at a point 
slightly displaced from its centre. To simplify the calculation, we place the dipole at the 
centre of the particle, and introduce a virtual propulsion force F p to effect propulsion. 
The magnitude Fo of this force is proportional to Sq. The directions of S p and F p are 
determined by p. In isolation, a self-propelled particle moves at a constant velocity u^p, 
with the speed uq determined by Sq. When it coexists with many such particles, its 
hydrodynamic interaction with the other particles alters its velocity and, more impor- 
tantly, its orientation. As a result, the motion of the particle is chaotic. Our simulations 
are not restricted to low particle concentration, as we implement the full hydrodynamic 
interactions between the particles, but we restrict the motion of particles to two dimen- 
sions to reduce computation. We have studied the statistical properties of a suspension 
of self-propelled particles for a range of the particle concentration, quantified by the area 
fraction </> a . We find several interesting features in the microstructure and statistics. We 
find that particles tend to swim in clusters wherein they are in close proximity. Con- 
sequently, incorporating the finite size of the particles and the near-field hydrodynamic 
interactions is of the essence. There is a continuous process of breakage and formation 
of the clusters. We find that the distribution of particle velocity at low and high 4> a are 
qualitatively different; it is close to the normal distribution at high a , in agreement with 
the experimental measurements of Wu fc Libchaberl (|2000l ). The motion of the particles 



is diffusive at long time, and the self-diffusivity decreases with increasing a . The pair 
correlation function shows a large anisotropic buildup near contact, which decays rapidly 
with separation. There is also an anisotropic orientation correlation near contact, which 
decays more slowly with separation. 



1. Introduction 

Nature presents a wide and fascinating array of organisms that can propel themselves 
in a fluid medium. Collections of swimming organisms exhibit intricate patterns and 
complex dynamics, such as the flocking of birds, schooling of fish and coherent motion 
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in microorganisms (jChildress et al. 19751 : Kessler 1986 : Wager 1911 ). While higher or- 
ganisms, such as fish and birds, have advanced sensory abilities to guide their motion, 
the sensory ability of microorganisms is quite rudimentary. Consequently, the interaction 
of microorganisms is largely mediated by the intervening fluid. Hence, understanding 
the hydrodynamics associated with the motion of individual organisms, and their fluid- 
mediated interactions is necessary for understanding their collective behaviour. 

The size a and swimming velocity Uq of most swi mming microorganisms are such that 
the R eynolds number Re = pu a a/rj is very small (|Lighthil]||l976l : |Purcel]||l977t iTavloi 
1951), and the Peclet number Pe = (Qnr]u^a) / '(hsT) is very large (jPedlev fc Kessler 
1992J). Here p and r\ are the density and viscosity, respectively, of the fluid, is the 
Boltzmann constant and T the absolute temperature. This means that the fluid motion 
is in the creeping flow regime, governed by Stokes equations, and Brownian motion of 
the microorganisms is negligible. If their density p p does not differ very much from p, as 
is normally the case, the Stokes number St = (p p / p)Re is also very small. In this regime, 
the inertia of the fluid and the 'particles' (i.e. the microorganisms) play no role, and 
hence propulsion does not come from bursts of acceleration generated by 'pushing' the 
fluid back, as in larger organisms. In Stokes flow, the net force on e ach swimmer is zero at 
every instant, and therefore the propulsion force balances the drag (Lighthill ll976HTavlorl 



195 lh . Instead, propulsion is achieved by a cyclic deformation of the body of the organism. 



The reversibility of Stokes equations implies that a reciprocal deformation during a cycle 
achieves no net displacement; hence a non-reciprocating, cyclic deformation is required. 

Microorganisms propel themselves in a number of ways: undulation of one or more 
flagella, helical motion of flagella, and coor dinated wavin g of a large number of cilia are 
some examples. Since the pioneering work of lTavlorl (|195ll ). a large number of st udies have 
consi dered the mechanics of propulsion by flagella and cilia (see, for example, Childress] 
1981), and a reasonable understanding of the subject has emerged. 

In this paper we focus on the collective behaviour of self-propelled particles in the 
regime of Stokes flow. For this purpose, we argue that the details of the mechanism of 
propulsion is not very important; regardless of the propulsion device, the fluid flow far 
from the particle is, to the lowest order of approximation, that due to a force dipole. 
We show that computing the hydrodynamic interactions between the self-propelled par- 
ticles is similar to that in a suspension of 'passive', or non-swimming, particles dis- 
persed in a fluid. It is well known that the hydrodynamic interaction between the 
suspended particles plays a crucial role in determining the bulk properties, such as 
its rheology, of a suspension. Moreover, the many-body hydrodynamic interactions re- 
sult in a range of complex behaviour , such as shear- induced diffusion and migration of 
particles dLeightqn fe Acrivodll987alfa: iNott fc Bradvlll994). anisotripi c microstructure 
( Bradv fe Morridll997l: IParsi fe Gadala-Marial|l987t ISingh fe Nottll2000h and non- linear 
rheology (|Sierou fe Bradvll2002l : ISingh fe NottJl200a IZarraga e£ aLll200oh . It is therefore 
quite likely that even our simple model will result in interesting and complex dynamical 
behaviour. 

Following the early work of Childress et al. ( 19751 ). several studies have considered the 
collective motion and pattern formation of populatio ns of self-propel led particles in a 



fluid by following a continuum mechanical approach. IChildress et ali showed that the 



'bioconvection' pa tterns observ ed in suspensions of motile organisms for over a century 
(see, for example, Wager|[l91l[ ) can be explained as a hydrodynamic instability akin to 
the Rayleigh-Benard instability. It is caused when the equilibrium between the negative 
geotaxis (i.e. their tendency to swim against gr avity) of the pa rticl es and their sedimen - 
tation due to their higher density is perturbed. Kesslerl d 1 986h andlPedlev et al\ (Il988[) . 
followed by other studies of the same group (jHill et al\ Il989t IPedlev fe Kesslerl Il990h , 
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coupled to this model the effect of gyrotaxis, or the competing effects of gravitational 
and viscous torques on the particles which together determine the swimming direction. 
In these continuum models, the system is modelled as a multiphase medium for which 
the field variables are the velocity and pressure of the suspension (fluid and particles), 
the number density of the particles and their orientation. The governing equations are 
the conservation of mass and momentum of the suspension, the number density of the 
particles and an evolution equation for the orientation. 



Mor e recently, another class of models has emerged, starting from the work of lVicsek et al 



(Il995l) . They proposed a model in which the position of the self-propelled particles evolve 
according to a simple set of rules: each particle moves with constant speed, and its orien- 
tation at any time step is the average orientation of other particles in its neighbourhood 
in the previous time step, with a random noise added. This simple model leads to a range 
of behaviour, including a continuous transition from a disordered 'phase' to an oriented 
phase with increasing number density and/o r decre asing noise. The continuum analogue 
of this model was presented by_ Toner fe Tul <l 19951). and i n a form more appropriate for 
freely swimming particles bv lSimha fe Ramaswamvl ( 20021) . the latter being an extension 
of the hydrodynamic theory of nematic liquid crystals. In all these models, there is an 
implicit assumption of the existence of short range forces between particles that result in 
alignment. These theories predict certain long wavelength instabilities and anomalously 
large fluctuations in the number density, which are yet to be tested experimentally. 

Though continuum models arc useful in understanding behaviour on large length and 
time scales, they cannot answer questions on the microstructure of the constituent par- 
ticles. Besides providing information on the small scale organization of the particles, 
knowledge of the microstructure also p rovides inputs to th e conti nuum models, and helps 
in refining them. The seminal work of iBatchelor fc Green 41972ft related the viscosity of 
a dilute suspension to the pair correlation of the suspended spheres; more recent studies 
dBradv fc Morrdll997t IParsi fc Qadala-Marialll987t ISierou fc Bradvll2002l ISingh fc NotH 
2000) have related the anisotropy of the particle microstructure to the non-linear rheology 
of the suspension, which is an important input in continuum models. 

In this study, we attempt to understand the collective motion of self-propelled par- 
ticles at a microscopic and mesoscopic level. We propose a model for self-propulsion, 
and incorporate the full hydrodynamic intera ctions interactions between the particles . 
We extend the Stokes ian Dynamics technique (jBradv et al. l ll988tlBradv fc BossdllQS^ 



Durlofskv et aZ.lll987l ) to incorporate self-propulsion, and carry out dynamic simulations 



for a range of the particle concentration. We track the motion of every particle as a func- 
tion of time, and extract the relevant statistical and microstructural properties. We find 
interesting and unexpected aspects of the collective dynamics reflected in the distribution 
of particle velocity and the position and orientation correlations. 

During the course of our investigation, a few papers have appeared in which mod- 



els for self-propulsion have been proposed ( Hernandez- Ortiz et al . 2005; Ishik awa et i_„ 
l2006t IRamachandran et al. 2006 ). Hernandez-Ortiz et al\ ( 20051 ) modelled a swimmer as 
a dumbbell comprising two beads connected by a rigid rod, with propulsion effected by 
a 'phantom' flagellum attached to one of the spher es. The force e x erted by the rigid 
rod is such that the net force on each bead vanishes. Ishikawa et al. ( 20061 ) developed a 
model of a squirmer, on the basis of the work of lLighthilll (|1952h . in which propulsion 
is generated by the tangen tial motion of the particle surface in a prescribed manner. 
IRamachandran et al. ( 20061 ) used the lattice-Boltzmann method (LBM) to simulate a 
swimmer, and achieved propulsion by an asymmetric distribution of forces on the sur- 
face of the particle with zero mean. In all these models, the net external force on the 
swimmer is zero, but there is a force dipole on it, as in our study. However, our model 
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differs from them in some significant ways: we consider the swimmers to be of finit e size 
and implement the full h ydrodynamic in teractions, while iHernandez-Ortiz et al\ treat 
the beads as point forces. Ilshikawa et al\ consider finite sized particles, but the nature 
of the model makes the computation of the hydrodynamic interactions between particles 
more difficult. However, it may be a more accu rate representation o f certain types of 
swimming organisms. While the LBM model of iRamachandran et al. also considers fi- 
nite sized particles, the computation of near-field interactions for particles in proximity is 
computationally intensive in their method (see the following paragraph) and inaccurate. 
Thus, we believe that the method we have proposed makes an optimal balance between 
accurate representation of the physical phenomenon on the one hand, and computational 
efficacy on the other. Moreover, the method we propose can be systematically refined by 
retaining higher moments of the surface force distribution, yielding a more sophisticated 
m odel for propulsion. 

IHernandez-Ortiz et al. ( 2005 ) used their model of a swimmer to simulate the collective 
dynamics of a system of particles bounded by plane parallel walls. Due to the nature of 
their model, they computed onl y the far field, point-forces int eractions between the beads . 
iLlopis fc Pagonabarragal (|2006l ) used the propulsion model of lRamachandran et all ([2006) 
to simulate interacting swimmers; however, they choose a particle radius of only 2.5 
times the lattice spacing, making the spatial discretisation very coarse. Moreover, they 
assumed an elastic collision between particles that are close to each other, which differs 
qualitatively from the dissipative lubrication interaction that is in force near contact. 
Our simulations show that a typical swimmer comes in close proximity with its neigh- 
bours; hence, we believe it is quite important to account for the finite particle size in 
the far field interactions, and accurately represent the near-field interactions to correctly 
capture the dynamics of the particles. We conside r the collective dynamic s of swimmers 
in a unbounded (spatially periodic) domain, while iHernandez- Ortiz et al\ studied a sys- 
tem bounded by plane parallel walls. Finally we present results on the microstructure 
and statistics of self-propelled particles that, to our knowledge, have not been reported 
earlier. 



2. Our model for a self-propelled particle 

To motivate our model for a swimmer, consider the schematic of a bi-flagellate or- 
ganism, shown in Fig. [TJ We ignore, for the moment, the effect of gravity or any other 
external force. The periodic, but non-reciprocating 'beating' of its flagella results in the 
action of a forward propulsive force ^F p by the fluid on each flagella. The motion of the 
particle is retarded by the fluid with drag force F p , as the particle cannot accelerate in 
the regime of Stokes flow. Though there is no net force, it is clear that there is a force 
dipole S p , and perhaps higher multipoles, acting on the particles. If there is no external 
torque on the particles, the dipole must be symmetric, i.e. it is a stresslet. We have used 
the biflagcllate, in which one can make the distinction between the 'propulsion arm' and 
the 'body' of the swimmer, only as an evocative example. Often it is not possible to 
separate the propulsion force and the drag, an example being a waving filament or sheet 
(|Tavlorlll95lh . Thus, it is more accurate to say that, at lowest order, a Stokesian swimmer 
is propelled by a force dipole. 

While the magnitude of the stresslet changes over the duration of a cycle, we are 
interested in the behaviour over time scales much larger than the period of a beat, and 
therefore assume the magnitude to be constant. However, the principal directions of the 
stresslet may vary in time, as interactions with other particles will cause the particle 
to rotate. The simplest model for propulsion, therefore, is a stresslet S p of constant 
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Figure 1. Schematic diagram of a biflagellate microorganism. The propulsion force on the 
flagella are matched by the drag force on the body of the organism. 



magnitude acting on the particles (Fig. [2|); indeed, this is the first approximation of a 
swimmer, regardless of the actual mechanism of its propulsion. 

The contribution of the stresslet on a swimmer to the bulk stress in the suspension has 
been recognised for long ( Pedlev fc Kessler 1990h . but it is not widely recognised that, 
at lowest order, the stresslet also that generates propulsion. 

The diameter of the flagella or cilia is typically far smaller than the body of the 
organism. For instance, the body diameter and flagella length of Chlamydomonas Nivalis, 
a bi-fiagellate alg a, are roughly 10/im, but the diameter of the flagella is only 0.1/xm 
(jMelkonianl Il992l) . Propulsion is generated because the flagella beat rapidly, so that 
their charac t eristic speed is much larger than that of the entire organism (~ 100 fim/s). 
Jones et al. (11994 used the resistance coefficients for a model flagellum provided by 
LighthillT i 19761 ) . and estimated that the body moves roughly a tenth of its diameter for 
each beat of the flagellum. Therefore, for the slow movement of the entire organism the 
hydrodynamic resistance of the flagella is only a small fraction of its total resistance, and 
to a first approximation may be neglected. 

We make the additional simplification that the particles are rigid spheres, as it sim- 
plifies the analysis and significantly eases computation. A stresslet acting at the centre 
of a sphere does not lead to movement, hence it must be displaced from the centre. It 
is clear from Fig. Q] that the centre of the dipole is not coincident with the hydrody- 
namic centre of the swimmer. Though the particles are treated as spheres, they possess 
an orientation which determines the direction of propulsion. Considering particle a, if 
p" is the unit vector id entifying its orientation, the propulsion stresslet acting on it is 
(|Pedlev fc Kesslerlll990l) 



S (p a p a -l/3 6), 



(2.1) 



where So is its magnitude, and S is the unit tensor. is traceless, as the trace con- 
tributes to the isotropic pressure of the fluid, which is arbitrary in an incompressible 
fluid. However, the induced stresslet S" (see can have a finite trac e; it arises from the 
interactions between particles and is related to the particle pressure (|Jeffrev et aLlll993 : 
Nott &: Bradvlll994l ). We note that Sq can be positive or negative: it is positive when 
the propelling arms pull the particle from the 'front', and negative when they push it at 
the 'rear'. Both cases occur in nature; Chlamydomonas is an example of the former, and 
spermatozoa an example of the latter. 
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Figure 2. Our model for a self-propelled particle. Here p is the unit vector identifying the ori- 
entation of the particle. Propulsion is generated by the action of a stresslet Sp = So (pp — 1/3 S) 
at a point displaced from the centre in the direction of p. 

From the linearity of Stokes equations, it follows that the velocity of locomotion of 
particle a is related to its propulsion stresslet Sp* as 

u^u^ + MZ-S;, (2.2) 

~ aa 

where M v s is the mobility of particle a due to the stresslet acting on itself, the caret 
denoting that it is the mobility for the off-centre stresslet, and u a is the velocity of the 

*• aa 

imposed macroscopic flow field at the particle centre. For an isolated particle, M us is a 
constant; in the presence of other particles (or boundaries), it depends on their positions 

~ aa 

relative to a. However, determination of M us is not a simple task: one way of doing 
it is to transfer the off-centre stresslet to the centre of the particle, which results in the 
introduction of all the higher multipoles at the centre. The dipole and all odd multipoles 
acting at the centre of the sphere do not lead to locomotion, as a result of symmetry, 
and hence one has to take account of the quadrupole and higher even moments. This is 
a level of detail we wish to avoid, as we would like to restrict the description to the level 
of monopole (force) and dipole (torque, and stresslet). To avoid this complication, we 

— aa 

resort to an artifice that makes the determination of M us unnecessary: we determine 
the propulsion velocity of particle a as though a virtual propulsion force F™ acts on it. 
This is an approximation, but we believe it to be a reasonable one; it only affects the 
way the propulsion velocity of a particle is hindered in the presence of other particles. 
Consequently, (|2.2[) is replaced by 

u* = u a + MZ a F -F^ (2.3) 

where My a F is the mobility of a sphere due to a force acting at its centre. As the force 
must act in the direction of p a , we set 

F°=F p a , (2.4) 

where F is the magnitude of the force. Clearly, F is determined by |So| (only the 
absolute value is relevant, as the direction of locomotion is the same whether the particle 
is being pushed or pulled), the particle radius a and the displacement of the stresslet 
from the centre; we therefore set 

F Q = \\So\/ ai (2.5) 
where A is an 0(1) dimensionless parameter that is related to the displacement of 
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from the centre. IPedlev fc Kesslerl (|l99(lh arrived at a relation similar to ()2.5() between 
the thrust exerted by the flagella an d the magnitude of the stressl et. In this simplified 
form, our model is similar to that of Hernandez-Ortiz et al. ( 2005f) . who used a 'phan- 



tom flagellum' to drive the dumbbells (see JQ). The important difference between our 
model and theirs is that interactions between particles modulate the effect of the virtual 
propulsion force, as the mobility Muf (see below) of a particle is reduced if it is close 
to another particle or a wall. 

We emphasise that the other particles do not perceive a force acting on a, but only 
the stresslet S^. The force on each particle serves only to determine its propulsion. 

In addition to the propulsion stresslet, a stresslet S" is induced on a when it is placed 
in a non-uniform velocity field, as a result of its rigidity. 

To summarise, our model for a self-propelled particle is the following: each particle a 
is treated as a rigid sphere with an orientation vector p a , on which a stresslet (of the 
form given in eq. !2.1[) acts at a point slightly displaced from its centre. The displacement 
of the stresslet from the centre of a is not perceptible to any other particle, and hence a 
appears as a sphere with a stresslet acting at its centre. The propulsion velocity of a is 
determined by introducing a virtual propulsion force F™ (of the form is given in eqs 12.41 
and !2.5|) on it. However, other particles do not perceive the force on a. 



3. Collective dynamics 

To place the problem in the framework of the Stokesian Dynamics method, let us first 
consider the dynamics of a suspension of passive particles. (We emphasise that the word 
'passive' is used in this paper to refer to particles that are not self-propelled, and not for 
tiny tracer particles that translate with the local fluid velocity.) For the motion of non- 
Brownian, rigid spheres in a Newtonian fluid at small Stokes number, Newton's second 
law reduces to 

F H + F = 0, (3.1) 

where F H is the vector of the hydrodynamic forces and torques on all the particles, and 
F is the vector of the non-hydrodynamic forces (external forces such as gravity, and 
inter-particle forces) and torques. 

To determine the motion of the particles, we must relate F K to the vector of their 
velocities and angular velocities it; in the creeping flow regime, this is accomplished by 
solving the Stokes equations, with the boundary conditions of no penetration and no 
slip on the surface of the particles. The hydrodynamic interactions between the particles 
cause the velocity and angular velocity of a given particle to depend not just on the 
hydrodynamic force and torque acting on itself, but also on the forces and torques acting 
on all the other particles. More precisely, u depends on the distribution of the hydro- 
dynamic force on the surfaces of all the particles. A convenient way of representing the 
distribution of force on the surface of a particle is the multipole moment expansion (see, 
for example, Kim &: Karrilal 1991 . chaps. 2-4): the zeroth multipole (monopole) is the 



net force on the particle; the first multipole (dipole) may separated into its symmetric 
and antisymmetric parts, the former being the torque and the latter the stresslet. By the 
linearity of Stokes equations, u is given by 

(3.2, 

where we have used (|3.1[i to replace F H with F. Here, S is the vector of the hydrodynamic 
stresslets, and u and e are vectors of the velocities and strain rates, respectively, of the 
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externally imposed flow field at the particle centres. The quantity Ai is the so-called 
'grand' mobility tensor, which can be decomposed into the mobility tensors representing 
the various couplings, 

the subscripts indicating the nature of the coupling. In principle, the right hand side 
of (|3.2|) must include higher moments of the force distribution on the particle surface, 
and the left hand side the higher gradients of the imposed velocity field — for each higher 
moment included, t here will be an additional equation. The Stok esian Dynamics method 
( Brady et a/]|l988t iBradv fc Bossislll988t iDurlofskv eTai\\l98l\ ), which we shall modify 



and use for the present problem, retains only the monopole and the dipole moments. 

The physical meaning of (|3.2p is as follows: the first line enforces (|3.1|1 . and determines 
u; the second line can be thought of as the equation to determine S. The stresslet is 
induced on each particle by the flow around it (externally imposed, and generated by the 
motion of other particles) so as to keep it rigid. 

The main advantage of the Stokesian Dynamics method is that Al (or equivalently 
the grand resistance 1Z = MT ) is computed and assembled in an accurate and efficient 
manner. This is done as a matched sum of the far-field and near-field interactions, 

M- 1 =M a 1 +Ti nf -n^. (3.4) 

Here, Adg is the mobility tensor that captures the far-field interactions, TZ D f is the 
near-field resistance tensor, and 72.^ is the far-field part of TZ D f that is subtracted to 
get a uniform as ymptotic expans i on. T hough Alff is assembled pair-wise, it has been 



demonstrated by IDurlofskv et al. ( 19871 ) that its inversion captures the many body hy- 



drodynamic interactions. 

The above framework must now be modified to incorporate the model for a swimmer 
we have developed in Jj2j In principle, the propulsion of the particles does not depend on 
the external force; hence, they can swim even when _F=0. However, we have introduced 
the virtual propulsion force F p in S}2] to avoid the inclusion of moments higher than the 
dipole in the multipole expansion. But F p must be recognised is a special force, as F£ 
acting on particle a only determines its velocity, and has no effect whatsoever on the 
other particles. The other particles perceive only the stresslet acting on particle a. In 
addition, the propulsion stresslet S p must be treated separately from the induced stresslet 
Si — the former is an inherent property of the swimmers, whose magnitude is constant 
and directions are determine by their orientations, whereas the latter is induced by the 
flow of the fluid around them. Accordingly, we distinguish the virtual propulsion force 
F p from the sum of external and inter-particle forces F ext , and the propulsion stresslet 
S p from the induced stresslet Si. The vectors F p and S p depend on the orientations of 
the particles through (|2.1|) and (|2.4| . 

After incorporation of the above modifications, the first line of (|3.2p takes the form, 

u - u = Mff-Fp + M UF -F cxt + M US -(S P + Si) (3.5) 

where Myp is the self mobility, i.e. the mobility of each particle due to a force on itself. 
For the particle pair a-(3, the self mobility is 

M s $ af3 = Mf F 5 a p, (3.6) 

where 8 af3 is the Kronecker delta. The first term on the right hand side of (I3.5[) provides 
propulsion; only the self mobility acts on F p , so that the virtual propulsion force on a 
particle has no effect on the other particles, as per our prescription. The second term 
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gives the velocity due to external and inter-particle forces, if any, and the third term is 
the velocity caused by the induced and propulsion stresslets. 
The modified form of the second line of (13.21) is 



- e = M EF -F cxt + Mes-S; + Ml°g- sclt -S p . (3.7) 

The first term on the right hand side gives the strain rate of the particles (relative to the 
externally imposed strain rate) caused by the external and inter-particle forces, and the 
second and third terms give the strain rate due to the induced and propulsion stresslets, 
respectively; since the particles are rigid, the three terms sum to — e. The physical in- 
terpretation of this equation is that each particle is imbedded in a flow field generated 
by the forces and stresslets on all the other particles, in addition to the macroscopic 
external flow, which induces a stresslet on it due to its rigidity. This equation determines 
Si, which when substituted in (|3 . 5[) yields the particle velocity vector it. 

Equations (|3.5[) and (|3.7p fully determine the collective dynamics of a system of self- 
propelled particles. Once it is determined for a particular configuration of the particles, 
their position and orientation are updated by time integrating 

dx 

- = «, (3.8) 

over a small time step At; the process is repeated until the desired duration of the 
simulation is reached. Here x is the vector of position and orientation coordinates of 
all the particles. The simulation is started with an initial configuration Xq; in all the 
simulations, the initial position and orientation of the particles were randomly assigned 
to achieve a uniform distribution. 

We consider the motion of self-propelled particles in the absence of any externally 
imposed flow, i.e. u = 0, e = 0. The particles were neutrally buoyant, so there was 
no net gravitational force on them. An inter-particle repulsive force of very short range 
was applied between particle pairs, as in previous studies using the Stokesian Dynamics 
method, in order to prevent overlap during the finite tim e-step integration of (|3.8[) . The 



form, strength and range of this force was the same as in Nott Sz. Brady ( 1994 ) 



It is convenient to scale all the variables in the following manner: stresslets by |5o|, 
forces by |5o|/a, distances by a, velocities by uq = |5o|/(67r?7a 2 ) and time by q/uq. The 
only adjustable parameter then is A; in all our calculations, we take A = 1. 



4. Results and Discussion 

Simulations of self-propelled particle suspensions were performed with periodic bound- 
ary conditions imposed in all directions to achieve an unbounded system. For a system of 
N particles in three dimensions, the velocity vector u is of dimension 6N, hence solution 
of the linear equation (|3.5p requires 0([6A] 3 ) computations at each time step. Similarly, 
(|3.7p requires 0([5A] 3 ) computations for the determination of the 5 N vector S[. As we 
require the simulations to run for a long duration to gather the statistical data of in- 
terest (see below), the computation requirement is considerable. To keep computation 
at a manageable level, we performed quasi-two-dimensional simulations, in which the 
particles were restricted to move in the (x,y) plane. For a fixed particle concentration, 
this reduces the numbe r of particles by a po wer of 2/3, and the sizes of u and Si to 3 AT. 
.Previous studies (e.g., iNott fc Bradvfli994) have shown the results of 2-d simulations 



to be similar to that of 3-d simulations, if the area fraction of the former is mapped 
suitably to the volume fraction of the latter. Nevertheless, it is desirable to study the 
motion in three dimensions of a large number of interacting particles, and we intend to 
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do so in a futur e inves tigation by using the Accelerated Stokesian Dynamics scheme of 



Sierou fc Bradvl (|200lh 



Most of our simulations were performed with 20 particles in a square unit cell of size 
L. The particle concentration, quantified by the area fraction </> a = Ntt/L 2 , was varied 
by changing L. A few simulations were performed with 30 and 40 particles (keeping a 
constant) to assess the effect of system size on the results. It was found that the effect 
of N on the statistical properties was quite small for </> a = 0.025, and imperceptible at 
higher </> a . Each simulation was run for 6000 dimensionless time units, but the data of 
the first 1000 time units was discarded for the statistical analysis, so as to ensure that 
the results are for a statistically steady state. 

Movies of simulations for </> a = 0.05 and 0.1 accompany this paper as supplementary 
material. Several features of the dynamical behaviour can be observed in the movies. (1) 
The orientation and velocity of the particles become randomised within a short period 
of time, and the motion of each particle resembles Brownian motion. (2) Though the 
particles are initially distributed uniformly (randomly) in the unit cell, in the dynamical 
steady state pairs, triplets and larger clusters are evident, within which particles are in 
close proximity. These clusters remain intact only for a short period of time; there is a 
continuous process of breaking up and formation of clusters. At any one instant of time, 
there are several groups of particles, a few stragglers, and relatively large empty spaces. 
(3) There is a tendency for particles that come in close proximity to align and move in 
such a way that one trails the other. (4) There is a substantial range in the velocity of the 
particles, from much lower to much higher than the free swimming velocity of a particle. 

Some of the features described above may be discerned from the snapshots shown in 
Fig. [3] As stated earlier, the initial configuration for this simulation was that of randomly 
assigned position and orientation of the particles. Each snapshot shows several pairs or 
larger clusters wherein the particles in close proximity. A striking example of a set of 
particles with like orientation moving in a train is seen in the last snapshot. 

We now proceed to analyse these features in greater detail. The results presented in 
for jj4.H -i j4.3l are for the case of 'pullers', i.e. So > 0. The case of 'pushers', i.e. So < 0, is 
discussed in MA\ 



4.1. Self diffusion 

The chaotic motion of the particles seen in the movies is a result of the hydrodynamic 
interactions between the particles. The interactions result in the perturbation of their 
velocity from their swimming velocity; more importantly, the vorticity generated by the 
motion of all the other particles causes each particle to rotate, thereby altering its orien- 
tation and hence its swimming velocity. The chaotic motion of the particles is apparent 
in Fig. |4l which shows the trajectories of three few particles in a particular simulation. 

Chaotic particle motion leads to diffusive behaviour at long time scales. The plot of 
the mean square displacement £ 2 = ((x(t) — Xo) 2 ) with time (Fig. [5]), the angle brackets 
indicating an average over many particles and initial conditions, shows that the motion 
is ballistic at small time (slope = 2), and diffusive at long time (slope=l). The time 
scale for transition from ballistic to diffusive motion decreases with increasing </> a , as per 
expectation, since interactions become stronger as </> a increases. We have determined the 
self-diffusivity from the Einstein relation 

£ 2 

V = lim — , (4.1) 

t^oo 4t 

and find that it is a decreasing function of the area fraction (Fig. [6|) . The inset of the 
figure shows that the diffusivity appears to obey a power law, V ~ <f>& n - This is in 
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Figure 3. Snapshots of particle position and orientation at various times in a simulation with 
particle concentration a = 0.05. The arrows indicate the orientation. Some particles appear to 
be in contact, but there is a thin lubricating layer of fluid between them. 
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Figure 4. The trajectories of three representative self-propelled particles in a simulation with 

particle concentration 4>a = 0.3. 
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Figure 5. The mean square displacement of the particles as a function of time. The dotted lines 
of slope 1 and 2 at the top are given to indicate the regimes of ballistic and diffusive motion. 



agreement with the results if IHernandez-Ortiz et al. ( 2005 ) , though the physical sig- 
nifi cance of a power lay decay is not clear to us. However, the diffusivities reported 
by IHernandez-Ortiz et al\ are much higher, except at very small For instance, at 



0.025 our diffusivity is just 20% lower than that reported bv lHernandez-Ortiz et al 



but at a = 0.025 it is lower by a factor of 10. Though the difference may be attributed 
partly to the fact that particle motion is restricte d to 2 dimensions in our s imula tions, the 
differences in the models must surely play a role: Hernandez-Ortiz et al. ( 2005[) imposed 
only the far field point-particle interactions, whereas our simulations have included the 
effect of finite particle size in the far-field interactions, and the strong near field inter- 
actions. The latter slow down particles that are in close proximity (but not similarly 
aligned) ; the pair distribution in H4.3I shows that there is a high probability of finding a 
p article in close proximity with another. 

[Wu fc Libchaberl (2000) conducted imaging experiments of Escherichia coli swimming 
in a freely suspended horizontal film, with a trace amount of passive spheres added. As 
the motion of the bacteria was restricted to the plane of the film, their experiments are 
similar to our simulations. However, they report the diffusivity of the passive spheres, 
and not the bacteria themselves. They find the diffusivity to be much larger than the 
Brownian diffusivity of the tracers, and using the Stokes-Einstein relation, they extract 
an "effective temperature" that is about 100 times greater than the temperature of the 
film. However, the effective temperature is not a useful quantity, as the diffusivity is only 
indirectly related to the temperature of the fluid. The temperature determines the rate at 
which the organism provides energy for locomotion, and the viscosity of the fluid, which 
together determine the swimming speed. It is more appropriate to scale the diffusivity 
by UQa, as we have done here. Scaled in this manner, the diffusivity they report for a 
4.5 y«m particle in a suspension of 1.8% by volume bacteria is w 1.6, while we find it to 
be about 0.8 for the swimmers. 

Self diffusion is also observed in sheared non-Bro wnian suspensions of passive spheres 
( Leighton fc Acrivosl Il987a ; Morris fc Brady 19961 ). where again velocity fluctuations 
arise from hydrodynamic interactions. However, the generation of fluctuational motion 
by changes in the particle orientation is not present there. 




4>a. 

Figure 6. The self diffusivity of the particles as a function of the area fraction. 



4.2. Distribution of particle velocity 

To study the distribution of particle velocity, one usually examines the probability density 
function f(u x , u y ), which is defined so that nf(u x , u y )du x du y is the probability of finding 
a particle whose x and y components of the velocity lie in the range du x and du y (around 
at u x and u y ), respectively. Here, n is the number density of the particles. We find it 
more convenient to study the velocity distribution in a particular direction, say the x 
direction, and define f x (u x ) as the probability density of the x velocity, regardless of the 
value of u y . The two functions are related by 



fx{u x ) = f(u x ,Uy)du y . (4.2) 



Let us first consider the case of a system of particles that are so far apart that they do 
not interact, i.e. a — > 0. In this situation, each particle moves with a constant dimen- 
sionless speed of unity in the direction of its orientation vector p. If the orientation of the 
particles is randomly distributed, the velocity distribution f m (u x ,u y ) (the superscript 
'NT indicating that the particles are non-interacting) is the Dirac delta function on the 
circle u r — 1, where u r = (u x + Uy) 1 / 2 is the speed , i.e. 

f m (u x ,u y ) = ±-6(u r -l). (4.3) 

The factor l/(2-7r) ensures that the integral of f m {u Xl u y ) over all velocities is unity. The 
distribution of the x velocity then follows from (|4.2[) , 

oo oo 

f x \u x ) = ^ J5(u r -l)du y = ^ JS(u r -l)duy. (4.4) 

-oo 

The latter equality arises from the symmetry of the integrand about u y — 0. As u x is 
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kept constant in the integral, du y — (u r /u y )du ri and hence 



C N1, \ 1 fS(u r -l)u r 



/*>*) = - jj^m^r. (4.5) 

From the sifting property of the delta function, we therefore have 

/fK) = I J, \, 1/2 (4-6) 

Note that f x l diverges as u x — > ±1. As the orientations are uniformly distributed, and 
there is no external force favouring motion in a particular direction, this distribution 
holds for all directions in the plane of motion. The distribution given by (|4.6j) is shown 
by the dotted line in Fig. [7J 

The distribution will, of course, be altered when the particles interact. We have anal- 
ysed the results of our simulations to determine the velocity distribution. The particle 
velocities were collected at dimensionless time intervals of unity, during which time a 
freely swimming particle moves a distance of its radius. The distribution function f x (u x ) 
was determined by constructing a histogram of the number distribution in equally sized 
intervals of u x , and normalising it so that J f x du x = 1. The same was repeated for the y 
direction. Due to the absence of directionality in the problem, the distribution function 
should be the same for all directions, which is indeed what we observe in Fig. [7j (compare 
the the lines and symbols). Given the isotropy of the velocity distribution, we henceforth 
denote by f(u) the distribution in any direction. It is evident from Fig.[jjthat the velocity 
distribution deviates from that of non-interacting swimmers (dotted line) even at small 

0a- 

For <fi a = 0.025, the velocity distribution is close to / NI for small but departs from 
it when |w| is greater than a value slightly less than unity — it has maxima at u sa ±1, and 
decays rapidly for \u\ > 1. Thus, there is a finite probability of finding a particle with a 
velocity significantly higher than that of an isolated swimmer. For large </> a , the velocity 
distribution is very different from / NI ; it appears to resemble the normal distribution 
(Fig.©, 



/(MH vib exp 



(«-«) 2 i (47) 



2a 2 



where u is the mean velocity and er 2 the variance. Small deviations from the normal 
distribution are apparent, such as a slight deficit of the prob ability at small M , and a 
faster decay at large |w| (see inset of Fig. [8]). In the study of Wu fc Libchaber ( 2000f ). 



referred to earlier, it is reported that the distribution of the speed u r of Escherichia coli 
follows the Maxwell distribution f(u r ) = exp[— u 2 /(2tr 2 )] at large particle concentra- 
tions, which is in accord with the normal distribution of the velocity in any direction. 
Since the largest concentration they studied is a volume fraction of = 0.1, it appears 
that they found a normal distribution at lower concentrations than we do. Whether the 
difference may be ascribed to the differences in the conditions of the experiments and the 
simulations, or the simplicity of our model is difficult to say. Our results indicate that 
careful measurements of the velocity distribution for a range of concentration is neces- 
sary. Conducting simulations and experiments in which particles are free to move in all 
three spatial dimensions would also be a worthwhile pursuit. Nevertheless, the qualitative 
agreement is perhaps an indication that our description is fundamentally sound, and it 
captures some of the important features of collective motion. 

Figure [9] compares the results obtained with and without the inclusion of the near- field 
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Figure 7. The probability distribution function of particle velocity in a suspension of self-pro- 
pelled particles at concentrations a = 0.025 and 0.3. The lines are symbols are the distributions 
of u x and u y , respectively, where x and y are the coordinates of the (fixed) laboratory refer- 
ence frame. The equality of f x and f y shows the absence of directionality in the problem. The 
dotted line is the velocity distribution for a collection of non-interacting, randomly oriented 
self-propelled particles, given by (|4.6|l . 
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Figure 8. The velocity distribution for a suspension with <f> a — 0.3 compared with the normal 
distribution having the same mean and variance. The inset shows the same plot in semi-log 
coordinates. 



hydrodynamic interactions, represented by the term TZ n f — in (|3.4| . We note that 
the velocity distribution is unchanged by the inclusion of near-field interactions for a 
dilute suspension of swimmers, but it is significantly altered for a relatively concentrated 
suspension. This is not an unexpected result, as the frequency with which a typical par- 
ticle comes into close proximity with others is relatively low at small </> a , but it increases 
with <f> a . It is pertinent to note that our simulations without the n ear-field interactions do 
not re duce to point-particle simulations, of the kind performed bv lHernandez-Ortiz et al. 
(l2005h . In Stokesian Dynamics simulations, the finite size of a particles is accounted for 
by retaining the induced dipole moments, and parts of the quadrupole and octupole mo- 
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Figure 9. Comparison of the velocity distribution for a dilute (0 a = 0.025) and concentrated 
(<^>a = 0.3) suspension of swimmers with and without the inclusion of near-field hydrodynamic 
interactions (NF). 
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Figure 10. The velocity distribution for a range of the particle concentration. 



ments (called the irreducible moments), in the multipole expansion of the force density 
distribution on the particle surface, and also the corresponding finite-size term s in the 



distribution on t he particle surrace, and also the c orresponding hmtc-size terrr . 
Faxen relations ( Brady et ofl ll988llBrad.v fc Bossidll98a iDurlofskv et aZJIl987h 



Reverting to our simulations with the full hydrodynamic interactions, the velocity 
distribution for all the particle concentrations that we have studied are shown in Fig.[T0l 
As a increases, the depth of the well between the two maxima decreases, vanishes 
completely at </> a slightly over 0.1, and the distribution resembles a normal distribution 
at high (j) & . The variance a of the distribution decreases with increasing </> a . 

4.3. Correlations 

As discussed earlier, the movies (see supplementary material) and the snapshots in Fig. [3] 
show significant correlation in the position and orientation of the particles. We first anal- 
yse the correlation in particle position in terms of the pair correlation function g{r2,T\), 
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Figure 11. A pair of neighbouring self-propelled particles. The position and orientation corre- 
lation functions are determined in a reference frame whose origin is coincident with the centre 
of particle 1, and whose x axis is in the direction of p 1 . 



which is defined so that ng(ri, r^) df*2 is the probability of finding particle 2 within the 
volume dr2 if particle 1 is situated at r\. As the system is spatially homogeneous, g is 
a function only of the separation r = r% — r%. As a result, in two dimensions we may 
express it as g(r, 9), where r is the scalar separation, and 9 an angle. It is not useful to 
measure 9 from a fixed laboratory axis, as the absence of directionality implies that g is 
isotropic. However, there is no isotropy if 9 is measured from the orientation vector of 
particle 1. This definition of 9 is also useful, as it tells us at what angle with respect to 
the orientation vector of a particle is there a greater likelihood of finding another. We 
therefore define 9 as the angle measured clockwise from p 1 to r, as shown in Fig [TT1 
Symmetry of the particle shape about its orientation axis results in the same symmetry 
for g, and hence g(r, 9) = g(r, 2tt — 9). We therefore consider the variation of g only for 
the first two quadrants, < 9 < it. 

Figure [T2l shows a greyscale plot of g(r,9). The strong anisotropic accumulation of 
particles near contact is apparent. There is a higher probability of finding another particle 
near its front (0 < 9 < tt/2) than near its rear (ir/2 < 9 < tt). The sharp decay of the pair 
correlation with separation distance is also apparent in the figure. This becomes clearer 
if we consider its angle averaged value g(r) = (l/7r) / g(r,9) d9, shown in Fig [13] Note 
that the pair probability near contact is far higher than that of a hard-sphere fluid at 
thermodynamic equilibrium. (While contact of smooth spheres is forbidden in Stokes flow, 
particles do come quite close to each other. For the purpose of this discussion, we do not 
distinguish between contact (r = 2) and near contact (2 < r < 2. 025)). For 4> a = 0.025, 
for ex ample, g(2) is just a little over unity for a hard-sphere fluid ([Carnahan fe Starling 
I1969T) . but here it is about 30 times larger. There is a similar difference in the build-up near 
contact for all particle concentrations. Secondly, g(r) decays much mor e rapidly with r 
than for a hard-sphere fluid, or a sheared suspension of passive particles ( Sierou fc Bradvl 
l2002h . Thus, the probability of finding a neighbour in close proximity is high, but it decays 
to the bulk probability within a short separation. 

It is pertinent to note that a very large build-up of particles near contact is seen 
in sheared suspensions of passive particles in the compressional quadrant, and the ac- 
tual value is found to be sensitive to the strength and range of the inter-particle re- 
pulsive force flS ierou fc Bradvl l2002HSingh fc Nottl [2000) , or the 'thermodynamic' force 
(jMorris fc Katvalll2002l ; iPhung et a.Z.Ml996f) that arises from Brownian motion. We have 
not varied the form of the repulsive force in this study, but believe that the results are 
relatively insensitive to it. The reason is that the balance between the hydrodynamic 
and repulsive (or thermodynamic) forces that exist in the compressional quadrant in 
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Figure 12. Greyscale plot of g(r,9) for cj> a = 0.05. The solid circle represents the surface of 
particle 1, r = 1, and the dashed circle the locus of centres of particle 2 if it were in contact 
with particle 1, r = 2. The radial distance beyond r — 2 has been stretched by a factor of 10, in 
order to discern the variation near contact. The arrow indicates the direction of the orientation 
vector p 1 . The bar below the plot gives the relation between the grey level and g. 
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Figure 13. The angle averaged pair correlation function as a function of the scalar separation, 
for three values of <j> a . The inset shows the variation at small r. 



sheared suspensions (|Bradv k Morris! 11997) is absent here; as two self-propelled parti- 
cles approach each other, they continue to rotate, and the difference in their orientation 
causes them to move apart. This process is clearly observable in the movies. As a result, 
we do not observe the long-lasting clusters that are seen in sheared suspensions. 

The angular variation of the pair correlation is determined by averaging over annular 
shells of width Ar = 0.025. The angular variation in the first and second shells are 
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Figure 14. Angular variation of the pair correlation function, averaged over annular shells of 
width Ar = 0.025. The upper cluster of lines are for the first shell 2 < r < 2.025, and the lower 
cluster of lines for the second shell 2.025 < r < 2.05. 

shown in Fig. [14] The plots for all the concentrations are qualitatively similar; they show 
a higher probability of finding a neighbour towards the front of each particle than in its 
rear, but only in the first shell. In the second shell, there is a smaller peak near 9 = 'in/ A 
in all the cases, which reflects the 'peeling off' of the accumulation from near contact, as 
seen in Fig. [TSl apart from this peak, the pair correlation at all angles is roughly uniform. 
In the third and higher shells, the pair correlation at all angles is roughly uniform, and 
close to the bulk value of unity (not shown). Thus, there is anisotropy in the distribution 
of the neighbours, but only at short separations. 

Finally, we consider the correlation in the orientation of particle pairs. For particles 1 
and 2 located at n and r 2 , respectively, the orientation correlation function is (p 1 -p 2 ), 
the angle brackets indicating an average over many particles and over time. It too is 
a function of r and 9, as defined in Fig [IT] A greyscale plot of (p 1 -p 2 ) is shown in 
Fig. [TSJ note that, unlike in Fig.[J2] the radial distance is not stretched here. The positive 
correlation of the orientations at the rear of the particle, around 9 = n, and negative 
correlation around 9 = 3ir/4 are evident. There appears to be no correlation at the front 
of the particle. 

To consider the variation of (p 1 -p 2 ) with r, we show its value as a function of r for 
9 = 7T and 37r/4 in Fig. [16] The positive correlation for the former and negative correlation 
for the latter near contact is evident. The correlation vanishes at large r, as expected, 
but its decay with r is much slower compared to that of g(r) (compare Fig. [T3|). meaning 
that particle orientations remain correlated over longer distances. The angular variation 
of (Pi'P 2 ) (Fig. fT7|) . averaged in an annular shell of width Ar = 0.1, shows a monotonic 
rise with 9, with maximum correlation near 9 = tt. 

Considering Figs. [T3l and IT4l we see that while the probability of finding a neighbour is 
high near the front of a particle, the probability of the neighbour being of like alignment 
is highest at the rear. 

4.4. Statistical features of a suspension of pushers (So < 0) 

As stated earlier, the results presented in §4.H - f4~3l are for So > 0, corresponding to case 
of the propelling arms of the swimmer pulling it from the front (see Fig.[lJ. As mentioned 
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Figure 15. Greyscale plot of the orientation correlation (p 1 -p 2 ). The solid circle represents the 
surface of particle 1, r — 1, and the dashed circle the locus of centres of particle 2 if it were in 
contact with particle 1, r = 2. The arrow indicates the direction of the orientation vector pj. 
The bar below the plot gives the relation between the grey level and (p 1 -p 2 )- 
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Figure 16. The two particle orientation correlation function ($>t"P-i) as a function of the scalar 
separation for two values of 6. The upper cluster of lines, showing positive correlation, are for 
8 — tt, and the lower cluster of lines, showing positive correlation, are for 9 — 37r/4. 



in 33 the opposite case of the propelling arms pushing from the rear, which in our model 
is achieved by setting So < 0, is also observed in nature. The collective dynamics of a 
suspension of such particles is therefore also of interest. 

We have performed simulations for this case at a particle concentration of a = 0.05. 
The mean-square displacement, diffusivity, velocity distribution, and the radial variation 
the pair correlation g(r) are found to be virtually identical to that of pullers (shown 
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Figure 17. Angular variation of the two particle orientation correlation function, averaged 

over the annular shell 2 < r < 2.025. 



in figures l5l [51 ITD1 and [T51 respectively), and are therefore not presented. The angular- 
variation of the pair correlation and the orientation correlation are, however, different. 
A greyscale plot of the pair correlation is shown in Fig. [18J the differences with Fig. [12] 
are apparent. Here, there is a strong accumulation of neighbours closer to the front 
(0 < 9 < 7r/4), spread over a slightly longer radial distance, and a weaker accumulation 
that is roughly uniform over other angles. As in Fig.[T21 the peeling off of the accumulation 
from near contact at the rear is evident. 

The greyscale plot of the orientation correlation (p 1 , p 2 ) (Fig- [HI) is also significantly 
different from the corresponding plot for pullers (Fig. [T5|) . Here, the region of positive 
correlation is spread over a larger range of 9 at the rear of the particle, and the region of 
negative correlation is pushed towards the equator (9 = it/2"). There is a second region of 
positive correlation just front of the equator, which was weaker and spread around 9 = 
in Fig. US 

This brings us to the question of whether the observed correlations, and the differences 
between the two cases, can be explained by simple mechanistic arguments, suc h as in a 
suspension of passive particles (jBatchelor fc Green 19721 : Bradv fc Morris! 1997 ). We are 
unable to provide a simple explanation, for the reason that as a swimmer approaches 
another, it is rotated by the fluid vorticity generated by the others, which changes its 
swimming velocity. Even in a dilute suspension of swimmers, where one may assume 
interactions to be pair-wise, the problem of determining the microstructure is signifi- 
cantly more complicated than in a suspension of passive particles: the trajectories of 
two particles, initially far apart, depend on their initial orientations. For determining 
the statistical properties of interest, their trajectories must be determined for all initial 
orientations, and the appropriate quantities averaged. 



5. Summary and conclusion 

We have studied the collective dynamics of self-propelled particles in a Newtonian 
fluid by conducting Stokesian Dynamics simulations. We have modelled each swimmer 
as a sphere whose propulsion arises by the action of a stresslet S p at a point slightly 
displaced from its centre. The strength So of the stresslet is assumed to be constant, and 
the principal directions of S p , and therefore the direction of propulsion, are determined by 
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Figure 18. Greyscale plot of g(r,9) for swimmers that 'push' from the rear, i.e. So < (see 
eg. I2.1j l, This figure should be compared with Fig. [12] for particles that pull from the front. 
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Figure 19. Greyscale plot of (p 1 -p 2 ) for swimmers that 'push' from the rear, i.e. So < (see 
eg. 12.111 . This figure should be compared with Fig. [T5] for particles that pull from the front. 



the orientation vector p of the particle. Rather than calculate the mobility due to an off- 
centre stresslet, we have determined the propulsion velocity of each particle by employing 
the ansatz of a virtual propulsion force F p acting in the direction of p. However, the force 
on a given particle only goes to determine its propulsion velocity, and all other particles 
only perceive the stesslet acting on it. 
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The chaotic motion of the interacting particles yields diffusive motion at long times, 
and the self diffusivity decreases a s the particle concentration <ft a increases. This trend 
is in agreement with the results of Hernandez-Ortiz et all ( 2005f ). but their diffusivities 
are up to an order of magnitude higher, as they ignored the near field hydrodynamic 
interactions. From the results of our simulations, we have extracted some important 
statistical indicators of the dynamics and microstructure. At high <p a , we find that the 
distribution of particle velocity u in any given direction i s close to the normal dist ribution, 
which is in accord with the experimental observation of lWu fc Libchaberl (|2000h . At low 
4>a,, however, the velocity distribution is qualitatively different: it has a local minimum at 
u = 0, peaks near ±wo, where uq is the speed of a solitary swimmer, and decays rapidly 
for larger velocities. 

Our analysis of the correlation of positions and orientations of particle pairs shows 
strong correlation near contact. The pair correlation function shows a large buildup of 
particles near contact even at low <j) a , suggesting that even at low particle concentration, 
the effects of finite particle size and the strong lubrication interactions are important 
in determining the collective dynamics. However, the pair correlation function decays 
much more rapidly with separation that for a hard-sphere fluid or a sheared suspension 
of passive (non-swimming) particles. Its angular variation shows anisotropy in the distri- 
bution of neighbours near contact, with a greater probability near the front of the test 
particle than in the rear. This anisotropy too vanishes quite rapidly with separation. The 
orientation correlation function decays relatively slowly with separation, and its angular 
variation shows greater correlation at the rear of the test particle than in the front. Thus, 
while there is a greater probability of finding a close neighbour at the front, the prob- 
ability that the neighbour has like orientation is highest at the rear. This result tallies 
with our observation in the movies that particles that come in close proximity often leave 
with like alignment, one trailing the other. 

A comparison of the statistical properties of suspensions of 'pullers' and 'pushers' reveal 
interesting similarities and some differences. The mean-square displacement, diffusivity, 
velocity distribution, and the radial variation the correlations are virtually identical in 
the two cases, but there are differences in the angular variation of the correlations. A 
mechanistic explanation for the nature of the correlations eludes us, as the rotation 
of particles as they approach each other makes their dynamics unamenable to simple 
analysis. 

A few words on how our model may be improved are in order. As remarked earlier, 
we consider this to be a simple 'first cut' model, which captures the most important 
features of interacting self-propelled particles. One important improvement would be to 
compute the correct mobility for an off-centre stresslet on a sphere, which necessitates 
the inclusion of higher moments of the force distribution. As a starting point, it appears 
useful to include a quadrupole at the particle centre; while a dipole at the centre of a 
sphere does not result in propulsion, due to symmetry, a quadrupole does. Another useful 
extension would be to consider non-spherical swimmers, in order to simulate the motion 
of organisms such as E. coli, which are rod-like. A simple way of extending the current 
framework to study rodlike swimmers is by 'sticking' two or more spheres together to 
form a linear extended object, and using constrained dynamics to ensure that they move 
as a solid body. Lastly, we note that while no external or propulsive torque was imposed 
on the particles in this study, it is straightforward to impose either. An external torque 
arises, for example in gravitaxis when mass is asymmetrically d istributed abo ut the centre 
of the the particle, and its alignment differs from the vertical (Kessler 19861) ; an internal 
or propulsive torque causes the 'tumbling', or sudden change in orientation, of bacteria 
like E. coli. 
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We have benefited from discussions with Sriram Ramaswamy and Ganesh Subramanian 
during the course of this work. 
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